function [neg_likelihood,nlh_grad]=nlh_garchmidas(parameter,g1,K,r,X)

% this is negative likihood function

miu=parameter(1);
alpha=parameter(2);
beta=parameter(3);
theta=parameter(4);
w=parameter(5);
m=parameter(6);


[T,C]=size(r);


% tau series construction -----------------------------------------------
tau = (m^2) + (theta^2) * ( X * betapolyn(K,[(K-1):-1:1]',w) ); % function of theta and w 
% -----------------------------------------------------------------------
tau_zero=find(tau<=0);
L_zero=length(tau_zero);
if L_zero >0
    tau(tau_zero);
    parameter
    stop
else
end;

% g series construction --------------------------------------------------
g=zeros(T,1);
g(1)=g1;
for i=2:T
    g(i)=(1-alpha-beta)...
         + alpha * ( (r(i-1)-miu)^2 ) / tau(i-1)...
         + beta * g(i-1);
end;
% ------------------------------------------------------------------------


% negative likelihood function --------------------------------------
neg_likelihood= T/2*log(2*pi) + (1/2)*sum( ((r-miu).^2)./(tau.*g) ) + (1/2)*sum( log(tau.*g) );
% -------------------------------------------------------------------

dg_dmiu=zeros(T,1);
dg_dalpha=zeros(T,1);
dg_dbeta=zeros(T,1);
dg_dtheta=zeros(T,1);
dg_dw=zeros(T,1);
dg_dm=zeros(T,1);

k_vec=[(K-1):-1:1]';

dtau_dtheta = (2*theta) * ( X * betapolyn(K,[(K-1):-1:1]',w) );

Nf=ones(1,(K-1))*(1-k_vec/K).^(w-1);
mf=sum( ((1-k_vec/K).^(w-1)).*log(1-k_vec/K) )/Nf;

dtau_dw     = (theta^2) * X * ( betapolyn(K,k_vec,w).*log(1-k_vec/K) )...
              - (theta^2) * X * betapolyn(K,k_vec,w) * mf;
dtau_dm     = (2*m);

for i=2:T
    dg_dmiu(i)   = -2*alpha*(r(i-1)-miu)/tau(i-1) + beta*dg_dmiu(i-1);
    dg_dalpha(i) = -1 + ((r(i-1)-miu)^2)/tau(i-1) + beta*dg_dalpha(i-1);
    dg_dbeta(i)  = -1 + g(i-1) + beta*dg_dbeta(i-1);
    dg_dtheta(i) = -alpha*((r(i-1)-miu)^2)/(tau(i-1)^2)*dtau_dtheta(i-1) + beta*dg_dtheta(i-1);
    dg_dw(i)     = -alpha*((r(i-1)-miu)^2)/(tau(i-1)^2)*dtau_dw(i-1) + beta*dg_dw(i-1);
    dg_dm(i)     = -alpha*((r(i-1)-miu)^2)/(tau(i-1)^2)*dtau_dm + beta*dg_dm(i-1);
end;

dL_dmiu   =(1/2)*( sum( -2*(r-miu)./tau./g ) - sum( ((r-miu).^2)./tau./(g.^2).*dg_dmiu ) + sum( dg_dmiu./g ) );
dL_dalpha =(1/2)*( -sum( ((r-miu).^2)./tau./(g.^2).*dg_dalpha ) + sum( dg_dalpha./g ) );
dL_dbeta  =(1/2)*( -sum( ((r-miu).^2)./tau./(g.^2).*dg_dbeta ) + sum( dg_dbeta./g ) );
dL_dtheta =(1/2)*( -sum( ((r-miu).^2)./(tau.^2)./(g.^2).*(dtau_dtheta.*g + tau.*dg_dtheta) ) + sum( (dtau_dtheta.*g + tau.*dg_dtheta)./tau./g ) );
dL_dw     =(1/2)*( -sum( ((r-miu).^2)./(tau.^2)./(g.^2).*(dtau_dw.*g + tau.*dg_dw) ) + sum( (dtau_dw.*g + tau.*dg_dw)./tau./g ) );
dL_dm     =(1/2)*( -sum( ((r-miu).^2)./(tau.^2)./(g.^2).*(dtau_dm.*g + tau.*dg_dm) ) + sum( (dtau_dm.*g + tau.*dg_dm)./tau./g ) );

% gradient ---------------------------------------------------------
nlh_grad=[dL_dmiu;dL_dalpha;dL_dbeta;dL_dtheta;dL_dw;dL_dm];
% ------------------------------------------------------------------
